﻿using Common.FrontEnd;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Numerics;
using System.Text;
using Util.Props;

namespace FrontEnd.Transform
{
    /// <summary>
    /// Computes the Discrete Fourier Transform (FT) of an input sequence, using Fast Fourier Transform (FFT). Fourier
    /// Transform is the process of analyzing a signal into its frequency components. In speech, rather than analyzing the
    /// signal over its entire duration, we analyze one <b>window</b> of audio data. This window is the product of applying a
    /// sliding Hamming window to the signal. Moreover, since the amplitude is a lot more important than the phase for speech
    /// recognition, this class returns the power spectrum of that window of data instead of the complex spectrum. Each value
    /// in the returned spectrum represents the strength of that particular frequency for that window of data.
    /// <p/>
    /// By default, the number of FFT points is the closest power of 2 that is equal to or larger than the number of samples
    /// in the incoming window of data. The FFT points can also be set by the user with the property defined by {@link
    /// #PROP_NUMBER_FFT_POINTS}. The length of the returned power spectrum is the number of FFT points, divided by 2, plus
    /// 1. Since the input signal is real, the FFT is symmetric, and the information contained in the whole vector is already
    /// present in its first half.
    /// <p/>
    /// Note that each call to {@link #getData() getData} only returns the spectrum of one window of data. To display the
    /// spectrogram of the entire original audio, one has to collect all the spectra from all the windows generated from the
    /// original data. A spectrogram is a two dimensional representation of three dimensional information. The horizontal
    /// axis represents time. The vertical axis represents the frequency. If we slice the spectrogram at a given time, we get
    /// the spectrum computed as the short term Fourier transform of the signal windowed around that time stamp. The
    /// intensity of the spectrum for each time frame is given by the color in the graph, or by the darkness in a gray scale
    /// plot. The spectrogram can be thought of as a view from the top of a surface generated by concatenating the spectral
    /// vectors obtained from the windowed signal.
    /// <p/>
    /// For example, Figure 1 below shows the audio signal of the utterance "one three nine oh", and Figure 2 shows its
    /// spectrogram, produced by putting together all the spectra returned by this FFT. Frequency is on the vertical axis,
    /// and time is on the horizontal axis. The darkness of the shade represents the strength of that frequency at that point
    /// in time: <p> <br><img src="doc-files/139o.jpg"> <br><b>Figure 1: The audio signal of the utterance "one three nine
    /// oh".</b> <p> <br><img src="doc-files/139ospectrum.jpg"> <br><b>Figure 2: The spectrogram of the utterance "one three
    /// nine oh" in Figure 1.</b>
    /// </summary>
    public class DiscreteFourierTransform: BaseDataProcessor
    {
        /** The property for the number of points in the Fourier Transform. */
        [S4Integer(defaultValue = -1)]
        public static String PROP_NUMBER_FFT_POINTS = "numberFftPoints";

        /** The property for the invert transform. */
        [@S4Boolean(defaultValue = false)]
        public static String PROP_INVERT = "invert";

        private Boolean isNumberFftPointsSet;
        private int numberFftPoints;
        private int logBase2NumberFftPoints;
        private int numberDataPoints;
        private Boolean invert;

        private Complex[] weightFft;
        private Complex[] inputFrame;
        private Complex[] from;
        private Complex[] to;

        private Complex weightFftTimesFrom2;
        private Complex tempComplex;

        public DiscreteFourierTransform( int numberFftPoints, Boolean invert ) 
        {
            this.numberFftPoints = numberFftPoints;
            this.isNumberFftPointsSet = (numberFftPoints != -1);
            this.invert = invert;
        }

        public DiscreteFourierTransform() {
        }

       /*
       /// (non-Javadoc)
        *
       /// @see edu.cmu.sphinx.util.props.Configurable#newProperties(edu.cmu.sphinx.util.props.PropertySheet)
        */
        override
        public void newProperties(PropertySheet ps) 
        {
            base.newProperties(ps);
            numberFftPoints = ps.getInt(PROP_NUMBER_FFT_POINTS);
            isNumberFftPointsSet = (numberFftPoints != -1);
            invert = ps.getBoolean(PROP_INVERT);
        }


        /// <summary>
        /// @see Sphincs.frontend.DataProcessor#initialize(Sphincs.frontend.CommonConfig)
        /// </summary>
        override public void initialize() 
        {
            base.initialize();
            if (isNumberFftPointsSet) {
                initializeFFT();
            }
        }


        /** Initialize all the data structures necessary for computing FFT. */
        private void initializeFFT() {
            /**
            /// Number of points in the FFT. By default, the value is 512,
            /// which means that we compute 512 values around a circle in the
            /// complex plane. Complex conjugate pairs will yield the same
            /// power, therefore the power produced by indices 256 through
            /// 511 are symmetrical with the ones between 1 and 254. Therefore,
            /// we need only return values between 0 and 255.
             */
            computeLogBase2(numberFftPoints);
            createWeightFft(numberFftPoints, invert);
            initComplexArrays();
            weightFftTimesFrom2 = new Complex();
            tempComplex = new Complex();
        }


        /** Initialize all the Complex arrays that will be necessary for FFT. */
        private void initComplexArrays() 
        {

            inputFrame = new Complex[numberFftPoints];
            from = new Complex[numberFftPoints];
            to = new Complex[numberFftPoints];

            for (int i = 0; i < numberFftPoints; i++) {
                inputFrame[i] = new Complex();
                from[i] = new Complex();
                to[i] = new Complex();
            }
        }


        /// <summary>
        /// Process data, creating the power spectrum from an input frame.
        /// </summary>
        /// <param name="input">the input frame</param>
        /// <returns>a DoubleData that is the power spectrum of the input frame</returns>
        private DoubleData process(DoubleData input)
        {
            /// Create complex input sequence equivalent to the real
            /// input sequence.
            /// If the number of points is less than the window size,
            /// we incur in aliasing. If it's greater, we pad the input
            /// sequence with zeros.
            
            double[] _in = input.Values;

            if (numberFftPoints < _in.Length) 
            {
                int i = 0;
                for (; i < numberFftPoints; i++) 
                {
                    inputFrame[i] = new Complex(_in[i], 0.0f);
                }
                for (; i < _in.Length; i++) 
                {
                    tempComplex = new Complex(_in[i], 0.0f);
                    
                    inputFrame[i % numberFftPoints] = Complex.Add
                            (inputFrame[i % numberFftPoints], tempComplex);
                }
            } 
            else 
            {
                int i = 0;
                for (; i < _in.Length; i++) 
                {
                    inputFrame[i] = new Complex(_in[i], 0.0f);
                }
                for (; i < numberFftPoints; i++) 
                {
                    inputFrame[i] = new Complex(0.0d,0.0d);
                }
            }

            
            // Create output sequence.
            
            double[] outputSpectrum = new double[(numberFftPoints >> 1) + 1];

            // Start Fast Fourier Transform recursion
            
            recurseFft(inputFrame, outputSpectrum, numberFftPoints, invert);

            
            /// Return the power spectrum
            DoubleData output = new DoubleData
                    (outputSpectrum, input.getSampleRate(),
                            input.getFirstSampleNumber());

            return output;
        }


        /**
        /// Make sure the number of points in the FFT is a power of 2 by computing its log base 2 and checking for
        /// remainders.
         *
        /// @param numberFftPoints number of points in the FFT
        /// @throws java.lang.IllegalArgumentException
         *
         */
        private void computeLogBase2(int numberFftPoints)
        {
            this.logBase2NumberFftPoints = 0;
            for (int k = numberFftPoints; k > 1; k >>= 1, this.logBase2NumberFftPoints++) 
            {
                if (((k % 2) != 0) || (numberFftPoints < 0)) 
                {
                    throw new ArgumentException("Not a power of 2: " + numberFftPoints);
                }
            }
        }


        /**
        /// Initializes the <b>weightFft[]</b> vector. <p><b>weightFft[k] = w ^ k</b></p> where: <p><b>w = exp(-2/// PI/// i /
        /// N)</b></p> <p><b>i</b> is a complex number such that <b>i/// i = -1</b> and <b>N</b> is the number of points in
        /// the FFT. Since <b>w</b> is complex, this is the same as</p> <p><b>Re(weightFft[k]) = cos ( -2/// PI/// k /
        /// N)</b></p> <p><b>Im(weightFft[k]) = sin ( -2/// PI/// k / N)</b></p>
         *
        /// @param numberFftPoints number of points in the FFT
        /// @param invert          whether it's direct (false) or inverse (true) FFT
         */
        private void createWeightFft(int numberFftPoints, Boolean invert) 
        {
            /**
            /// weightFFT will have numberFftPoints/2 complex elements.
             */
            weightFft = new Complex[numberFftPoints >> 1];

            /**
            /// For the inverse FFT,
            /// w = 2* PI / numberFftPoints;
             */
            double w = -2* Math.PI / numberFftPoints;
            if (invert) {
                w = -w;
            }

            for (int k = 0; k < (numberFftPoints >> 1); k++) 
            {
                weightFft[k] = new Complex(Math.Cos(w* k), Math.Sin(w* k));
            }
        }

        /// <summary>
        /// Reads the next DoubleData object, which is a data frame from which we'll compute the power spectrum. Signal
        /// objects just pass through unmodified.
        /// </summary>
        /// <returns>the next available power spectrum DoubleData object, or null if no Spectrum object is available</returns>
        override public IData getData()
        {
            IData input = getPredecessor().getData();

            getTimer().start();

            if ((input != null) && (input is DoubleData)) 
            {
                DoubleData data = (DoubleData) input;
                if (!isNumberFftPointsSet) {
                    /*
                    /// If numberFftPoints is not set by the user,
                    /// figure out the numberFftPoints and initialize the
                    /// data structures appropriately.
                     */
                    if (numberDataPoints != data.Values.Length) {
                        numberDataPoints = data.Values.Length;
                        numberFftPoints = getNumberFftPoints(numberDataPoints);
                        initializeFFT();
                    }
                } else {
                    /*
                    /// Warn if the user-set numberFftPoints is not ideal.
                     */
                    if (numberDataPoints != data.Values.Length) {
                        numberDataPoints = data.Values.Length;
                        int idealFftPoints = getNumberFftPoints(numberDataPoints);
                        if (idealFftPoints != numberFftPoints) {
                            Trace.WriteLine("User set numberFftPoints (" +
                                    numberFftPoints + ") is not ideal (" +
                                    idealFftPoints + ')');
                        }
                    }
                }
                input = process(data);
            }

            // At this point - or in the call immediatelly preceding
            // this -, we should have created a cepstrum frame with
            // whatever data came last, even if we had less than
            // window size of data.

            getTimer().stop();

            return input;
        }


        /**
        /// Returns the ideal number of FFT points given the number of samples. The ideal number of FFT points is the closest
        /// power of 2 that is equal to or larger than the number of samples in the incoming window.
         *
        /// @param numberSamples the number of samples in the incoming window
        /// @return the closest power of 2 that is equal to or larger than the number of samples in the incoming window
         */
        private  static int getNumberFftPoints(int numberSamples) 
        {
            int fftPoints = 1;

            while (fftPoints < numberSamples) {
                fftPoints <<= 1;
                if (fftPoints < 1) {
                    throw new Exception("Invalid # of FFT points: " + fftPoints);
                }
            }
            return fftPoints;
        }


        /**
        /// Establish the recursion. The FFT computation will be computed by as a recursion. Each stage in the butterfly will
        /// be fully computed during recursion. In fact, we use the mechanism of recursion only because it's the simplest way
        /// of switching the "input" and "output" vectors. The output of a stage is the input to the next stage. The
        /// butterfly computes elements in place, but we still need to switch the vectors. We could copy it (not very
        /// efficient...) or, in C, switch the pointers. We can avoid the pointers by using recursion.
         *
        /// @param input           input sequence
        /// @param output          output sequence
        /// @param numberFftPoints number of points in the FFT
        /// @param invert          whether it's direct (false) or inverse (true) FFT
         */
        private void recurseFft(Complex[] input,
                                double[] output,
                                int numberFftPoints,
                                Boolean invert) 
        {

            double divisor;

            /**
            /// The direct and inverse FFT are essentially the same
            /// algorithm, except for two difference: a scaling factor of
            /// "numberFftPoints" and the signal of the exponent in the
            /// weightFft vectors, defined in the method
            /// <code>createWeightFft</code>.
             */

            if (!invert) {
                divisor = 1.0;
            } else {
                divisor = numberFftPoints;
            }

            /**
            /// Initialize the "from" and "to" variables.
             */
            for (int i = 0; i < numberFftPoints; i++) 
            {
                to[i] = new Complex(0.0d,0.0d);
                from[i] = new Complex(input[i].Real/divisor,input[i].Imaginary/divisor);
            }

            /**
            /// Repeat the recursion log2(numberFftPoints) times,
            /// i.e., we have log2(numberFftPoints) butterfly stages.
             */
            butterflyStage(from, to, numberFftPoints, numberFftPoints >> 1);

            /**
            /// Compute energy ("float") for each frequency point
            /// from the fft ("complex")
             */
            if ((this.logBase2NumberFftPoints & 1) == 0) 
            {
                for (int i = 0; i <= (numberFftPoints >> 1); i++) {
                    output[i] = squaredMagnitudeComplex(from[i]);
                }
            } 
            else 
            {
                for (int i = 0; i <= (numberFftPoints >> 1); i++) {
                    output[i] = squaredMagnitudeComplex(to[i]);
                }
            }
      }

        /**
      /// Method to compute the squared magnitude of a complex number.
       *
      /// @return the squared magnitude of the complex number
       */
        public double squaredMagnitudeComplex(Complex _this)
        {
            double squaredMag;
            squaredMag = _this.Real* _this.Real + _this.Imaginary* _this.Imaginary;
            return squaredMag;
        }

        /**
        /// Compute one stage in the FFT butterfly. The name "butterfly" appears because this method computes elements in
        /// pairs, and a flowgraph of the computation (output "0" comes from input "0" and "1" and output "1" comes from
        /// input "0" and "1") resembles a butterfly.
        /// <p/>
        /// We repeat <code>butterflyStage</code> for <b>log_2(numberFftPoints)</b> stages, by calling the recursion with the
        /// argument <code>currentDistance</code> divided by 2 at each call, and checking if it's still > 0.
         *
        /// @param from            the input sequence at each stage
        /// @param to              the output sequence
        /// @param numberFftPoints the total number of points
        /// @param currentDistance the "distance" between elements in the butterfly
         */
        private void butterflyStage(Complex[] from,
                                    Complex[] to,
                                    int numberFftPoints,
                                    int currentDistance) {
            int ndx1From;
            int ndx2From;
            int ndx1To;
            int ndx2To;
            int ndxWeightFft;

            if (currentDistance > 0) {

                int twiceCurrentDistance = 2* currentDistance;

                for (int s = 0; s < currentDistance; s++) {
                    ndx1From = s;
                    ndx2From = s + currentDistance;
                    ndx1To = s;
                    ndx2To = s + (numberFftPoints >> 1);
                    ndxWeightFft = 0;
                    while (ndxWeightFft < (numberFftPoints >> 1)) {
                        /**
                        /// <b>weightFftTimesFrom2 = weightFft[k]   </b>
                        /// <b>                      *from[ndx2From]</b>
                         */
                        weightFftTimesFrom2 = Complex.Multiply
                                (weightFft[ndxWeightFft], from[ndx2From]);
                        /**
                        /// <b>to[ndx1To] = from[ndx1From]       </b>
                        /// <b>             + weightFftTimesFrom2</b>
                         */
                        to[ndx1To] = Complex.Add
                                (from[ndx1From], weightFftTimesFrom2);
                        /**
                        /// <b>to[ndx2To] = from[ndx1From]       </b>
                        /// <b>             - weightFftTimesFrom2</b>
                         */
                        to[ndx2To] = Complex.Subtract
                                (from[ndx1From], weightFftTimesFrom2);
                        ndx1From += twiceCurrentDistance;
                        ndx2From += twiceCurrentDistance;
                        ndx1To += currentDistance;
                        ndx2To += currentDistance;
                        ndxWeightFft += currentDistance;
                    }
                }

                /**
                /// This call'd better be the last call in this block, so when
                /// it returns we go straight into the return line below.
                 *
                /// We switch the <i>to</i> and <i>from</i> variables,
                /// the total number of points remains the same, and
                /// the <i>currentDistance</i> is divided by 2.
                 */
                butterflyStage(to, from, numberFftPoints, (currentDistance >> 1));
            }
       }

    }
}
